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ABSTRACT 

The N-independence observed in the evolution of massive black hole binaries (MBHBs) in 
recent simulation of merging stellar bulges suggests a simple interpretation beyond complex 
time-dependent relaxation processes. We conjecture that the MBHB hardening rate is equiva¬ 
lent to that of a binary immersed in a field of unbound stars with density p and typical velocity 
tr, provided that p and tr are the stellar density and the velocity dispersion at the influence ra¬ 
dius of the MBHB. By comparing direct N-body simulations to an hybrid model based on 
3-body scattering experiments, we verify this hypothesis: when normalized to the stellar den¬ 
sity and velocity dispersion at the binary influence radius, the N-body MBHB hardening rate 
approximately matches that predicted by 3-body scatterings in the investigated cases. The ec¬ 
centricity evolution obtained with the two techniques is also in reasonable agreement. This 
result is particularly practical because it allows to estimate the lifetime of MBHBs forming in 
dry mergers based solely on the stellar density profile of the host galaxy. We briefly discuss 
some implications of our finding for the gravitational wave signal observable by pulsar timing 
arrays and for the expected population of MBHBs lurking in massive ellipticals. 

Key words: black hole physics - galaxies: kinematics and dynamics - galaxies: evolution - 
gravitational waves - methods: numerical 


1 INTRODUCTION 

Massive black holes (MBHs) are fundamental building blocks in 
the process of galaxy formation and evo lution; they are ubiqui - 
tous in nearby galaxy nuclei (see, e.g., iMagorrian et aljfiggsh . 
and their masses corre late with the properties of the host bulge 
dKormendv & Hdl2013l . and reference therein). If MBHs are com¬ 
mon in galaxy centres at all epochs, as implied by the fact that 

! :alaxies harbour active nuclei for a short period of their lifetime 
Haehnelt & Rees|[l993h . then a natural consequence of the hier¬ 
archical paradigm of structure formation is that a large number of 
massive black hole binaries (M BHBs) forms along th e cosmic his¬ 
tory, following galaxy mergers dBegelman et ahlflgSOh . 

As galaxies build-up, they keep turning cold gas into stars, 
with the result that massive galaxies at low redshift are pre¬ 
dominantly gas poor. Nonetheless, massive gas poor galaxies 
cont inue to merge with satellites in groups and clusters (see, 
e.g. [van Dokkund l2005t) . In fact, a s ignificant fraction of them 
is observed in close pairs (see, e.g. iLopez-Sanjuan et aH l2012l : 
IXu et aljl2012l) . implying a significant merger rate. The fate of 
the forming MBHB in these massive, gas poor galaxy merg- 


* E-mail: asesana@stai'.sr.bham.ac.uk 
t E-mail:khanfazeel.ist@gmail.com 


ers have attracted the attention of several inve stigators that t ack¬ 
led the problem either (s emi)analvticallv (e.g. lOuinlanlfl?^: 1^ 
I2002L ISesana et al.l 20^^ or by mean s of direct N-body sim- 
ulations d Milosavl ievic & Merritt! 1200 ll : iMakino & Funatol |2004 
iBerczik et alJl20od. e.g.). 


In recent years, advancement in massive parallel computing 
made possible to simulate ’ab initio’ the evolution of MBHBs in 
merg i ng bulges, using up to a few million particles dPreto et al.l 
|201 ll: iKhan erZll201 iL EoTl: lOualandris & Merrittll2012h . These 
authors found that the MBHB hardening rate (the rate at which its 
semimajor axis shrinks) is independent of the number of particles, 
suggesting efficient diffusion of stars in phase space, keeping the 
binary ’loss cone’ - i.e. the portion of stellar distribution phase 
space collecting orbits with small angular momentum, intersecting 
the binary semimajor axis - full. 


In this letter we shall assume that these simulations capture 
the relevant physics governing the diffusion of stars in phase space 
in the merger remnant, and we interpret this in light of simple 3- 
body scattering theory. By comparing direct N-body simulations 
to an hy brid model ba sed on 3-body scattering experiments (devel¬ 
oped bv ISesan al2010l) . we show that, when normalized to the stellar 
density and velocity dispersion at the binary influence radius, the 
N-body MBHB hardening rate matches that predicted by 3-body 
scattering. This result is extremely practical, because the MBHB 
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coalescence time can then be reliably estimated once the density 
profile of the stellar distribution of the host galaxy is known. 

The paper is organized as follows. In Section 2, we describe 
the N-body runs, the hybrid model and the match between the two. 
We present our results in Section 3 and we discuss its importance 
and possible applications in Section 4. 


2 MATCHING N-BODY SIMULATIONS TO THE 
HYBRID MODEL 


2.1 N-body runs 


The N-body runs performed f or this study, closely follow those 
described in jKhan et al.l 120121 hereinafter K12). The two merg- 
ing bulges are initialized as equal mass and size Dehnen profiles 
(lDehneiilll993h : 


p{r) 


(3 - 7)M« rp 

dvr rT'(r + ro)^ ’’’ 


( 1 ) 


where M* is the total mass of the stellar bulge, ro is the scale ra¬ 
dius, and 7 is the inner logarithmic slope. A MBH is placed at the 
centre of each bulge. We perform four simulations, by assuming 
two different inner cusp slopes - 7 = 1,1.5 - and MBHB mass 
ratios - q — M 2 /Mi = 1, 1/3, where Mi > M 2 are the masses 
of the primary and secondary MBH respectively. In all cases we 
set Ml = O.OOSMoai. The two bulges are initially at a separation 
15 on an bound orbit with eccentricity 0.75. The units adopted for 
the integration are G = M, = ro = 1. A-body simulations are 
performed using />—GRAPE-l-GPU on accre, a high-performance 
GPU computing cluster at V anderbilt University , Nashville, TN. 
The code is upd ated version oflHarfst et alj i2007h and is described 
in section 2.2 of I Khan eta L H2013h . 


2.2 Hybrid model 


The h ybrid model we use has been extensively described in dSesanal 
I 2 OIOI ' hereinafter SIO). For a given density profile - e.g. the one 
given in equation Q- one can compute the initial separation oo 
at which the mass enclosed in the binary is twice the mass of the 
secondary (i.e. M*(< ao) — 2 M 2 ). At this point, an initial ec¬ 
centricity eo is also assumed, and the binary is evolved forward in 
time. The evolution takes into account for the initial scattering of 


bound stars l eading to the erosio n of the central stellar cusp (see 
full details in ISesana et alj|200§l . followed by a phase dominated 
by s cattering of unbound stars intersec ting the binary semimajor 
axis jOuinlanll 19961 : [^sana et al.l |2006h. and the efficient g ravita- 


tional wave (GW) emission stage 1 Peters & Mathewill96^ lead¬ 
ing to final coalescence of the system. After the short phase of cusp 
disruption, in the unbound scattering phase, the binary hardening 
proceeds at a rate 


d 

dt 




( 2 ) 


where H^h ~ 15 — 20 is a dimensionless rate, and p and a are 
effectively free parameters that, in the original formulation of the 
3-body scattering problem, represent the density and velocity of 
the distribution of intruding stars at infinity. The fact that Hsh is 
found to be independent of a (for hard binaries) means that the 
hardening rate given by equation ID) is about constant in time. In the 
framework of the hybrid model, one can tune the supply of unbound 
stars to the binary by picking some specific value of p and cr related 


to the adopted stellar distribution - equation (D- Units are set so 
that M = Ml + M 2 = G = ao = 1 and all the technical details 
can be found in S 10 . 


2.3 Comparing the MBHB evolution in the two cases 


As shown in the previous subsections, the N-body and the hybrid 
approaches have been developed using different units, which makes 
a comparison between the two somewhat cumbersome. One can 
assume the density profile given by equation Q and then con¬ 
vert the units of one model into the other, but this is not the best 
way to proceed, because the initial density profile of the N-body 
model evolves during the merger, and the stellar distribution in 
which the binary forms and evolves is different from the initial 
one. In some sense, we want to normalize the binary evolution 
to the stellar environment. Because of the fairly constant harden¬ 
ing rate seen in N-body simulations, it is natural to also write it 
in the form given by equation where one has the freedom to 
pick an arbitrary value of cr and p, which results in a dif ferent value 
of j^Nb (this approach has alread y been employed bv iKhan et al.l 
I 2 OI itlOualandris & Merritj 2012 ll . 

We conjecture that the A-body MBHB hardening rate is 
equivalent to that of a binary immersed in a field of unbound stars 
with density pinf and velocity crinf, equal to the density and veloc¬ 
ity dispersion of the surrounding stellar distribution at the binary 
influence radius. The latter is defined as the radius ,rinf, containing 
twice the binary mass in stars - i.e., M* (< rinf) = 2M. This state¬ 
ment is equivalent to say that the binary loss cone is full at rinf (see 
discussion in Section 3.1 in SIO). The explanation of why the loss 
cone is full is beyond the scope of this letter, and it must trace back 
to efficient stellar diffusion in the time-dependent, triaxial, rotating 
potential of a merger remnant. If our conjecture is correct then by 
plugging the {d/dt){l/a) measured in the N-body simulation, pinf 
and CTinf into equation 0 , we should obtain TTnu ~ Hoh. This 
also means that the evolution of the MBHB in the two approaches 
should be remarkably similar if the same time units are used. Here, 
for ’same time units’ we mean units that normalize the stellar den¬ 
sity and velocity in equation 0 to the same value in the N-body run 
and in the hybrid model. In practice the N-body and hybrid codes 
have their own units in which: 


d 

dtoh 

(-) 

V“3b/ 

crab 

(3) 

d 1 

djsiht ' 

-) 

vUNb/ 

GpNb „ 

= -nNb- 

CTNb 

(4) 


Therefore, assuming as a reference the N-body time, the MBHB 
evolution obtained in the hybrid formalism will match that of the 
N-body runs if the time in the hybrid model is rescaled according 
to: 


fflNb O-Nb p3b , 

iNb =-tab, 

flSb crab pNb 


(5) 


where cr^b = crinf and pNb = pinf, and all quantities are measured 
in the respective units. 


3 RESULTS 

3.1 Hardening rates 

A comparison between the four N-body runs and the hybrid models 
featuring the same parameters is shown in figure [T] The evolution 
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Figure 1. MBHB evolution in the N-body runs and in the hybrid model normalized to the N-body units. T = 0 has been set to the moment when the binaiy 
reaches the separation a = 0. In each column of plots we show (from the top to the bottom) the time evolution of the semimajor axis, the inverse of the 
semimajor axis and the eccentricity. In each panel, the green curve is the result of the N-body simulation, the red curve is the evolution predicted by the hybrid 
model with the same parameters. 


is rescaled to N-body units, i.e. the time unit in the hybrid model 
has been converted according to equation i). Moreover the time 
has been shifted so that T — 0 when the MBHB sits at a semimajor 
axis oq. We recall that ao is defined so that M*(< oq) = 2 M 2 , 
and it is the starting point of the hybrid model. We see in the top 
panels that, in all cases, the evolution of the semimajor axis matches 
quite well in the two approaches, as expected. Differences between 
the N-body runs and the hybrid model become more visible when 
1/a is plotted (central panels). Here we see that we still have a 
good match, even though in the N-body runs a somewhat slower 
evolution is visible. A slight decline of the shrinking rate can be 
due to either (i) a not completely full loss cone or (ii) a declining 
density of the underlying stellar distribution at the influence radius. 

To gain insights on this point, we plot in figure [2] the stel¬ 
lar density profile of the B3 N-body run at different snapshot (we 
should keep in mind that in figure [T] we normalized T to the mo¬ 
ment when the binary becomes bound, which corresponds to snap¬ 
shot 0062 in figure[2l(. In run B3, rinf is located at 0.084 (see table 
[T](, and we notice that at this separation the stellar density declines 
constantly with time. The drop ranges between 10% and 40% de¬ 
pending on the simulation, naturally slowing down the MBHB evo¬ 
lution by an amount comparable to what seen in figure|T] 

The hardening rate comparison is shown in table [T] where 
Hj^h has been computed by fixing rinf at the time the MBHB first 
becomes bound (column 2 in the table) and averaging over the en¬ 
tire subsequent evolution in hard binary phase. As expected, we 
find Htih ~ Hsb (within 30%), confirming our conjecture. We no¬ 
tice, however, that H^ih is systematically lower, likely due to the 
evolution of the density profile discussed above. In fact, in comput¬ 
ing T^Nb we fixed rinf at T = 0. However the subsequent density 
decline (figure |3, also implies a sizeable expansion of rinf dur¬ 
ing the simulation. In practice rinf is an evolving quantity, and any 
computation of a time averaged N-body hardening rate is intrinsi¬ 
cally approximate. Therefore, fixing rint{T — 0) is expected to 



Figure 2. Density profile of simulation B3 at different times (see labels in 
figure). The cyan dash-dotted line is the density profile of the initial model 
( 7 i = 1.5 and the black double-dashed line is the best fit to the density 
profile when the binary becomes bound (7 = 1 . 1 ). 


Model 

7i 

1 

ri 

ao 

Hnu 

Hsh 

AO 

1.0 

1 

0.16 

0.113 

14.53 

17.4 

A3 

1.0 

1/3 

0.13 

0.066 

14.29 

16.4 

BO 

1.5 

1 

0.10 

0.066 

11.43 

16.5 

B3 

1.5 

1/3 

0.084 

0.039 

12.46 

17.7 


Table 1. Hardening rates for the 4 simulated MBHB systems (column 1). 
We report the initial inner slope of the density profile (column 2), the binary 
mass ratio (column 3), the binary influence radius and initial separation of 
the hybrid model in N-body units (columns 4 and 5), and the hardening 
rates estimated from the N-body simulations and from the hybrid model 
(columns 6 and 7). 
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return an average -f^Nb that is lower to the one found in 3-body 
scatterings. 

Another remarkable aspect of this comparison is the quite 
similar evolution of the binary eccentricity in the two approaches 
(lower panels in figure [T), even though the N-body data are in¬ 
evitably noisy. We caution, however, against some of the limita¬ 
tions of the hybrid model. One major caveat is that the bound cusp 
is unaffected by the presence of the binary until M 2 reaches aq. 
Although this is a fair assumption for binaries with q < 1/10, 
in the comparable mass cases investigated here M 2 will affect the 
distribution of stars around Mi earlier in its evolution, implying a 
smoother transition between dynamical friction and scattering of 
stars refilling the loss cone. As a consequence, the binary evolu¬ 
tion in the short initial cusp erosion phase in the hybrid model is 
in general faster than in the N-body simulation. This is not a major 
issue, since the binding energy transferred to the binary from the 
cusp erosion is always the same, whether this happens promptly 
or gradually. However, the angular momentum exchange can be 
significantly different, meaning that the eccentricity evolution pre¬ 
dicted by the hybrid model might not be fully trustworthy in this 
phase. 


3.2 Coalescence times 


Although we plan to address extensively this topic in a forthcoming 
paper, the results obtained in the previous section can be used to get 
an approximate estimate of typical lifetimes and characteristic sep¬ 
arations of MBHBs forming in dry galaxy mergers. The evolution 
of the semimajor axis can be written as 


where 


da 

dt 


da 


da 


B 


“77 +“77 ——Aa -g 

dt 136 dt lew a'^ 


O'inl 


B = 


64G^MiM2MF(e) 

5c® 


( 6 ) 

(7) 


7?eff « 0.75ro(2i+®"+ - 1)“^ (lDehnenlll993h . This latter is 
a function of the galaxy stellar mas s which depends on t he na- 
ture of the galaxy host. In particular IPabringhausen et^ ( 1200Sh 
found that Reft/pc = max(2.95M*'f®®, 34.8M*for ellip¬ 
tical galaxies, whereas Rgg/pc = 2.95M* |®® for the bulges of 
spirals and ultra compact dwarfs (M*_6 is the mass of the stellar 
bulge normalized to 10® Mq). 

With the relations described above, for each MBHB mass, we 
can estimate M*, tq, Tint, Pint, Cinf of the host galaxy, and compute 
its lifetime and characteristic separation through equations l[^ and 
(|9). Results are shown in figure [3] for a range of MBHB and host 
galaxy properties. One of the consequences of our findings is that 
MBHBs in purely stellar environment do not stall, but their life¬ 
times are rather long, ranging from 0.1 to several Gyr, depending 
on the properties of the systems (left panels in figure]^. In partic¬ 
ular, binaries with M = 10® — 10^® Mq hosted in giant ellipticals 
with shallow density profiles, which should dominate the GW sig¬ 
nal in th e nHz frequency re gime accessible to pulsar timing arrays 
(PTAs) dSesana et al.ll2008h . might have lifetimes as long as 2-10 
Gyr, regardless of their mass ratio. These lifetimes, however, can 
be significantly shorter if eccentricity is driven to very high values 
during the stellar dominated evolution phase (as some of the tracks 
shown in figure [T]might suggest). For example, for e = 0.99, typi¬ 
cal coalescence times are a factor of « 20 shorter (lower left panel 
in figure[3]l. Regardless of their mass ratio or of the galaxy density 
profile most of this systems spend the majority of their time at a 
separation which is roughly proportional to the MBHB mass, and 
for M = 10® — 1O^®M0 systems is in the range 0.1-1 pc. High 
eccentricity promotes efficient GW emission at larger separations, 
and binaries with e = 0.99 besides having a lifetime which is « 20 
shorter, have a better chance to be found at separations which are 
« 10 times larger (lower right panel in figure[3}. 


4 DISCUSSION AND CONCLUSION 


and F{e) = (1 - e ®)~+^fl + (73/24)e® -f (37/96)e'‘] 

dPeters & Mathewall963h . Since the stellar hardening is oc a/ and 
the GW hardening is oc binaries spend most of their time 
at the transition separation obtained by imposing {da/dt)^^ = 
(da/dt)^^-. 


/gw 


64GVinfMiM2MT’(e) 


1/5 


5c^Hpinl 


(8) 


and their lifetime can be written as 


-a/ \ *+int 


(9) 


In fact the lifetime estimated through equation © differs by only 
about 10% from that obtained by integrating equation lO. To get 
physical estimates for realistic galaxies we need to estimate Cinf 
and Pint. We get the former fro m the M — a relation as re¬ 
ported bv lKorrnendv & h 3 d2013h . Mg = O.SOOctIoo*. where Mg 
is the binary mass normalized to 10 ®Mq and crgoo is the bulge 
velocity dispersion normalized to 200 km s“^. The latter is ob¬ 
tained from the density profile given by equation (ID evaluated 
at the influence radius rinf = ro/{[M*/(2M)]^+® — 1}, 

where rg and M* need to be specified. For a given MBHB mass 
M, we get M« from the M — bulge relation as reported by 
iKormendv & Hoi d2013h . Mg = 0.49M/'if, where where M*,ii 
is the bulge mass normalized to 10^^ Mq. The break radius rg is 
connected to the bulge effective radius Reg through the relation 


We performed the first detailed comparison between direct N-body 
simulations of MBHB mergers in stellar bulges and the hybrid 
model based on 3-body scattering experiments developed by SIO. 
Guided by the fairly constant in time, N-independent behaviour of 
the MBHB hardening in the N-body runs, we conjectured that the 
N-body MBHB hardening rate is equivalent to that of a binary im¬ 
mersed in a field of unbound stars with density pinf and velocity 
cTinf, equal to the density and velocity dispersion of the surrounding 
stellar distribution at the binary influence radius. We demonstrated 
the validity of this statement by showing that the dimensionless N- 
body hardening rate is comparable to that found in standard 3-body 
experiments, if normalized to pinf and (Tinf (see table [T}. We also 
showed that, when normalized to the same reference density and 
velocity dispersion, the MBHB evolution in the N-body runs and 
in the hybrid model matches fairly well, even in terms of eccentric¬ 
ity growth. The N-independent MBHB hardening seen in N-body 
runs, although fairly robust, has been tested up to 10® particles only. 
IVasilievl ( l2015h developed a clever technique to extrapolate N-body 
results In the formal limit N 00, finding hardening rates that are 
about a half of those found at W = 10® (IVasilietll2014l) . How¬ 
ever, these models are equilibrium triaxial systems, that might not 
capture non-relaxation and rotation effect present in merger rem¬ 
nants. Moreover, efficient diffusion in realistic systems due to the 
abundant presence of massive perturbers (globular clusters, giant 
molecular clouds) ca n easily bring the hardenin g rate back to the 
full loss cone regime dPerets & Alexand^l2008h . 
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Binary mass [Solar masses] 


Figure 3. Lifetime (left panels) and characteristic separation (right pan¬ 
els) of MBHBs vs binary mass. In the top panels solid-black, dashed- 
blue and dotted-red lines are for density profiles with 7 = 0.5,1,1.5 re¬ 
spectively; in the central panels solid-black, dashed-blue and dotted-red 
lines ai'e for MBHBs with q = 1,0.32,0.1 respectively; in the bottom 
panels solid-black, dashed-blue and dotted-red lines are for MBHB with 
e = 0,0.9,0.99 respectively. In each bifurcated curve, the upper branch 
is for regular' ellipticals, the lower branch is for bulges of spirals and ultra 
compact dwarfs. 


Our findings have a variety of interesting applications and im¬ 
plications in several astrophysical contexts. Firstly, one can directly 
infer the lifetime of a putative MBHB residing in a gas poor galaxy 
solely based on the observed stellar density profile and on an es¬ 
timate of the MBHB mass. In fact, those are the only ingredients 
needed to work out the the stellar density at the binary influence 
radius and to estimate the binary coalescence time through the hy¬ 
brid model. This prescription can be reliably used in semianalytic 
models or large scale simulations of galaxy formation, where often 
the binary hardening phase is by-passed, assuming prompt coales¬ 
cence. Secondly, the coalescence timescales we found for typical 
elliptical galaxies are longer than 1 Gyr, unless binaries are driven 
to considerably large eccentricities (e > 0.9). Together with the 
fact that massive ellipticals had pro bably undergone at least one 
major merger since z = 1 (see, e.g.. lvan Dokkumll2005l) . this im¬ 
plies that parsec scales MBHBs should be quite common in massive 
gas poor galaxies and their observational signatures might be ac¬ 
cessible by future 30-40 meter optical facilities. Moreover, binary 
lifetimes in the Gyr range might lead to the freq uent formation of 
MBH triplets following two subsequent mergers l lHoffman & LoebI 
l2007h . which can have impo rtant consequences for the expected 
GW signal in the PTA band jAmaro-Seoane et'^l2010l) . and we 
plan to investigate this point in future work. Thirdly, our results 
highlight the importance of knowing the eccentricity of the MBHB 
at formation. After the binary becomes hard, the eccentricity evo¬ 
lution can be reliably tracked using the hybrid model, but the value 
of the eccentricity when the binary first becomes bound cannot be 
predicted on the basis of simple arguments and it might well be 
stochastic (in fact in the K12 runs there is a large spread in the 


values of initial MBHB eccentricity). Highly eccentric binaries can 
coalesce much faster than their circular counterparts, having im¬ 
portant implications on the MBHB population, on the occurrence 
of triplets, and on the nature of the GW signal in the PTA band, as 
discussed above. A large set of fairly short simulations following 
the MBHB only to right after it becomes bound would be extremely 
useful to complete the picture, and we are currently pursuing this. 
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